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ABSTRACT 

We construct a simple symplectic map to study the dynamics of eccentric orbits 
in non-spherical potentials. The map offers a dramatic improvement in speed over 
traditional integration methods, while accurately representing the qualitative details of 
the dynamics. We focus attention on planar, non-axisymmetric power-law potentials, 
in particular the logarithmic potential. We confirm the presence of resonant orbit 
families (“boxlets”) in this potential and uncover new dynamics such as the emergence 
of a stochastic web in nearly axisymmetric logarithmic potentials. The map can also 
be applied to triaxial, lopsided, non-power-law and rotating potentials. 


1. Introduction 


The morphology of orbits in triaxial potentials determines the structure of triaxial galaxies. 
The simplest triaxial potentials are those of Stackel form, which support four major families of 
regular orbits (boxes, short-axis tubes, and inner and outer long-axis tubes) and no stochastic 
orbits (Lynden-Bell 1962, de Zeeuw 1985). Many plausible triaxial potentials with smooth cores 
exhibit similar structure: most of phase space is occupied by these four major orbit families, 
with the leftover phase space hosting resonant minor families and stochastic orbits (Schwarzschild 
1979). 

Potentials without cores, such as triaxial logarithmic potentials of the form ln(X] /af), can 
have quite different behavior (Gerhard &: Binney 1985, Miralda-Escude &: Schwarzschild 1989). In 
this case the box orbits are replaced by a number of minor orbit families corresponding to various 
resonances. The minor orbits are centrophobic (center-avoiding) in contrast to the centrophilic 
box orbits, which all pass arbitrarily close to the center. 

In special cases non-axisymmetric power-law potentials can have Stackel form. Sridhar &: 
Touma (1997) constructed planar, non-axisymmetric cuspy potentials which are separable in 
parabolic coordinates. These potentials remain separable if a central black hole is added. In the 
absence of a black hole, the dynamics are scale-free and all orbits belong to the centrophobic 
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“banana” family. With a black hole, an additional family of centrophilic orbits (“lenses”) emerges 
to replace the box orbits of models with a smooth core. 

The relevance of power-law potentials is enhanced by recent observations. High-resolution 
Hubble Space Telescope photometry of nearby elliptical galaxies and spiral bulges shows that the 
luminosity density near the center is a power-law, ^(r) oc with —0.3 ^ a ^ 1.7 (e.g. Gebhardt 
et al. 1996). Few or no galaxies have smooth cores {a = 2). There is also growing evidence 
for massive black holes in the centers of many nearby galaxies (Kormendy & Richstone 1995), 
which generate singular potentials that share many properties with the potentials generated by 
power-law densities. 

We would like to understand how the orbit structure in non-axisymmetric power-law potentials 
oc changes with the exponent a, and with the degree of non-axisymmetry. An exhaustive 
study with conventional integration methods is costly and difficult to interpret. We show how 
the dynamics of eccentric orbits in power-law potentials can be illuminated with the help of a 
symplectic mapping. The mapping models the evolution of such orbits as a two-step process: (i) 
precession of the orientation of the orbit in an axisymmetric potential; (ii) a kick to the angular 
momentum of the orbits at apocentre (where the star spends most of its time, and the torques are 
likely to be strongest). We use this mapping to study the orbital structure of non-axisymmetric 
power-law potentials over the range of a relevant to galaxies. 


1.1. Scale-free spherical potentials 


We assemble some properties of orbits in attracting spherical potentials of the form 


$a(r) = 


iFr“, a 7 ^ 0, 
Klnr 0 = 0. 


( 1 ) 


where we assume henceforth that —1 < a <2, which is true for most plausible potentials (a = — 1 
is Keplerian; a = 2 is the harmonic oscillator; a = 0 is the singular isothermal sphere). Potentials 
with smaller a arise from density distributions with greater central concentration. We shall call 
potentials with a > 0 concave potentials and those with a < 0 convex. With no loss of generality 
we can set K = sgn(Q;) for a ^ 0 and K = 1 for a = 0. The energy corresponding to a circular 
orbit with radius Tc is then 

J i|a|r“-bsgn(a)r“, a ^ 0, 

( ^ Inrc, a = 0. 

Note that sgn(£') = sgn(a) for all non-zero a > —2, so that 


\E\ = (l + ia) r“. 


a 7^ 0. 


( 3 ) 



The angular momentum of a circular orbit is 


'a 2 


= /i(q;)|£'|« + 2 J a / 0, 


LrXE) = \ l“l^^^(l + 2«) 

exp(i?) = /i(0)exp(£') 


exp 


a = 0. 


( 4 ) 


Now consider motion in a plane and let L be the scalar angular momentum, positive for prograde 
orbits and negative for retrograde orbits; it is convenient to work with the dimensionless angular 
momentum 

( 5 ) 


Lc{E) ’ 

which can vary from —1 to +1. The orientation of an eccentric orbit can be specified by the 
azimuthal angle of its apocenter, (pn- We may write 


(pn+i — (pn T 5(0, y); 


( 6 ) 


where the precession rate g is independent of energy because the potential is scale-free, and 
|y(a,y)| is 27rPr/P0, where Pr and are the radial and azimuthal periods. The function g is odd 
in y, and is given by 


dr 


g{a,y) = 2h{a)y \ ^ r[2ssn{a){r^ - rp - y^h^a^^ 


r[—2r^ Inr — y^/i^(0)]^'^^ ’ 


, a / 0, 


a = 0, 


( 7 ) 


where the integral is over all radii for which the radicand is positive. For near-radial orbits we 
have 



lim, g(a, y) = 

±50 (a), 

(8) 


y- 


where 








for a > 0, 



go{oi) 

= < 27r 

1 

for a < 0. 

(9) 

When y 

is small we have (cf. Appendix A) 




5(a,y) - ffo(a) 

~ \y\^ 

(10) 

where 







' 1 

1 < a < 2 , 



/3(a) = < 

a 

0 < a < 1, 

(11) 


1 2a 

TTFf ■ 

-| < a < 0, 


1 —1 < a < — 


For near-circular orbits we have 
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which is the usual epicyclic approximation. There are two special cases for which the precession 
rate is simple, 

g{a, y) = sgn(y) | (13) 

in addition, forQ; = |,l,—i,—| the precession rate can be expressed in terms of elliptic functions 
(Whittaker 1959). In general g{a,y) cannot be determined in closed form, although an asymptotic 
series is given in Appendix A. A useful identity, easily derived from (^), is (Grant & Rosner 1994) 

£/(«, y) = 9{a, y) where a = (14) 

thus the behavior of g in the range — 1 < a < 0 is determined by its behavior in the range 

0 < a < 2. 

The map we discuss in the following section repeatedly employs the function g{a, y) at fixed a. 
To minimize calculations of the integral (|^ we tabulate log]^Q[( 7 (a, y) — go{a)] at logio(y) = —6(0.1)0 
and interpolate using a cubic spline. Computing this integral is straightforward but difficult 
to do well, in part because of the integrable singularities at the endpoints and in part because 
the endpoints themselves are determined implicitly; we can supply our program or results upon 
request. Figure s shows g{a,y) - go{a). 


2. A map for eccentric orbits 


We now examine the behaviour of orbits in a scale-free axisymmetric potential that is 
perturbed by a small non-axisymmetric potential. If the unperturbed orbit is nearly circular, 
the behaviour can be analyzed using the epicycle approximation (e.g. Binney & Tremaine 1987, 
§3.3.3). Here we focus instead on an approximation valid for eccentric orbits. 

In most cases the torques exerted by a non-axisymmetric potential on an eccentric orbit are 
concentrated near apocenter, since (i) the lever arm is larger; (ii) the particle spends most of its 
time near apocenter; (iii) any external tidal forces are stronger at larger radii. 

Let yn and 4>n be the values of the dimensionless angular momentum and the azimuthal 
angle of a particle at its apocenter passage. If the non-axisymmetric potential has exp(im(/>) 
symmetry then the time-integrated torque over the apocenter passage can be written 
—eLc{E) sinm(/>n; the minus sign is appropriate if the azimuthal minimum of the non-axisymmetric 
potential lies along the ray 0 = 0. This torque induces a change in the angular momentum, of 
which half occurs before apocenter and half after; thus the dimensionless angular momentum after 
the particle leaves apocenter is 


yn = yn - ^esinm0n- 


(15) 
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Fig. 1.— Precession rate g{a,y) (eq. ^ as a function of dimensionless angular momentum y, for 
power-law potentials with exponent a = —0.8(0.2)1. 8 . The precession rate is an odd function of 
y and is only shown for y > 0. Potentials with a < 0 are represented by dashed lines, those with 
a > 0 by solid lines. The value of <7 at y = 0 (eq. P) has been subtracted off; thus the analogous 
curves for a = —1 and 2 are zero everywhere (cf. eq. [l^). 
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The position of the following apocenter is 

4>n+i = (t>n +g{a,y'J, (16) 

and the angular momentum at this apocenter is 

Vn+I = y'n-lesmm(l)n+i. (17) 

Equations (|l^-p^ define a simple map {(f>n,yn) —> {4’n+i,yn+i) that describes the dynamics 
of eccentric orbits. The map is symplectic, like the exact equations of motion, and is more 
general than the exact equations because it does not depend on the specific radial form of the 
non-axisymmetric potential, so long as the torque is concentrated near apocenter. 


2.1. Relation to the standard map 


When 1 < a < 2 equations i)-® imply that g{y) ~ tt + Cy for |y| -C 1, where C is a 
constant (Fig. |T] shows how well this approximation holds). Thus the mapping (|T^)- (|T^) reduces 
to 


y'n = yn-\^smm(l)n, 

</’n+l = 4>n ~\~ Cy'^, 

yn +1 = y'n-^esinmcpn+i- (18) 

If in addition m is even, we can change variables to q = mcp + n, p = mCy, to get 

Pn = Pn + ^K sin Qn, 

Qn+1 — Qn ~i~ Pni 

Pn+l = p'n + \K sin qn+l, (19) 


where K = mCe. In this form, the map is the symmetric expression of the Chirikov-Taylor map, 
otherwise known as the standard map. This map serves dynamicists as a laboratory for examining 
the behavior of area-preserving maps as one increases the strength of perturbations. We refer 
the reader to the vast literature on the standard map (e.g. Lichtenberg &: Lieberman 1992), and 
simply recall the transition to global stochasticity that occurs as K exceeds the critical value of 
~ 1, or equivalently when 

1 


Thus we expect that the map is regular when the non-axisymmetric perturbation is small, but only 
in the range 1 < a < 2. More specifically, we shall find below that when a < 1 near-radial orbits 
can exhibit a complex chaotic structure, even for arbitrarily small non-axisymmetric perturbations. 
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2.2. Non-axisymmetric potentials 


We now ask what scale-free non-axisymmetric potentials can be generated by plausible density 
distributions. Scale-free density distributions can be written in the form 

p{r) = aemYemiO, </>)■ (21) 

i,m 

The corresponding potential is (e.g. Binney &: Tremaine 1987, §2.4) 

When £ = 0, the range of validity of equation (|2^ can be extended to all a > 0, since the dynamics 
are determined by the radial force —d^/dr, which can be determined from (p^ whenever — 1 < a. 
Outside the range in which equation (p^ is valid, the multipole potential is determined by the 
distribution of material at very large or very small radii, so that it satisfies Laplace’s equation and 
a = ^ or —i — 1. Thus the range of potentials in which we are interested is given by 


-l<a, i = 0, 

-l-i<a<e, i>0. 


(23) 


Since we restrict ourselves to the range of exponents — 1 < a < 2, these constraints are always 
satisfied for monopole, quadrupole and higher multipoles (£ = 0 or £ > 2). For dipole potentials 
(£ = 1), the constraints are only satisfied in the smaller range — 2 < a < 1. 

To assess the realism of the map, we shall compare trajectories in the map to orbits in 
scale-free non-axisymmetric potentials of the form 


^a(xi,X2) 


/ a;2\ 2" 

sgn(a) ( x| j , a/0, 

i log , a = 0, 


(24) 


where xi = rcos4> and X 2 = rsincj) are the usual Cartesian coordinates. The potential is specified 
by the exponent a and the axis ratio of the equipotentials, h < 1. The maximum angular 
momentum for an orbit of given energy now depends on the azimuthal angle 4>, and is given in 
terms of the maximum angular momentum in the analogous axisymmetric potential (eq. iby 


Lmax(L', </>) — 


Lc{E) 


(cos^ cf) -|- sin^ 


thus \y\ = 


< (cos^ 0-|-sin^ (/>/5^) (25) 


The map that approximates motion in this potential is specified by the same a, azimuthal 
wavenumber m = 2, and the torque amplitude e. To relate e to the axis ratio h = 1 — 5, we examine 







near-radial orbits {y <C 1) in nearly axisymmetric potentials (5 -C 1). Then it is straightforward 
to show that to lowest order in 6 and y 


^ — Qa — ' 




t)'\ 


(l -h 


1 -^/ 3 ) 


h+h r(i + i/a) 
r(| + l/a)' 

(27re)^/2^ 

-H r(-| + i//3) 


r(i//3) 


a > 0, 
a = 0, 

a = —(3 < 0. 


(26) 


This ratio varies smoothly from Q 2 = it = 3.1416, to Qo = 4.1327, to Q-i = 27r = 6.2832. 

The map is based on the approximation that torques are concentrated near apocenter—that 
the intervals of the orbit when torques are exerted (mostly near apocenter) and when the orbit 
processes (mostly near pericenter) are disjoint. This approximation is most accurate for nearly 
radial orbits—if the orbit is almost exactly radial then all of the precession occurs as the orbit 
passes the origin. If the non-axisymmetric component of the potential is given by (24), then the 
map is most accurate for concave potentials (a > 0), since in this case the torque is concentrated 
at large radii. 

A second possible form for the potential is 


4>„(xi,X2) = c{xl - xj) + 


sgn(a) (xf-bx|)2", a^O, 
I log {xj + xl) , a = 0; 


(27) 


this form is not scale-free, but may be more natural if the non-axisymmetric potential arises from 
external tidal forces. For this potential the torques are always concentrated at large radii. 


2.3. The logarithmic potential 


We hrst apply the map to the logarithmic potential (a = 0), which is relevant to realistic 
galaxy models, and also represents a boundary between potentials in which near-radial orbits 
process by vr in one orbit and those that process by smaller angles (eq. |^. The behaviour of orbits 
this potential has been examined by Richstone (1980, 1982), Binney &: Spergel (1982), perhard 


m 


fc Binney (1985)| , Pfenniger fc de Zeeuw (1989)1 , and Miralda-Escude Sz Schwarzschild (1989) 


Figure 1^ shows the surface of section (SOS) for orbits in the potential {‘M) with a = 0 and 
b = 0.88. The SOS shows the dimensionless angular momentum and azimuth each time the orbit 
passes through apocenter. The dashed lines show the maximum dimensionless angular momentum 
allowed by (^). The analogous map, with azimuthal wavenumber m = 2 and dimensionless 
strength e = Qo{l — h) = 0.5, is shown in Figure ^3. The principal differences are: (i) the shapes 
of the regular orbits and the chaotic zones are slightly different; (ii) small resonant islands for 
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Fig. 2.— (a) Surface of section (SOS) for orbits in the potential CM), with axis ratio b = 0.88. 
Each orbit is plotted for 1000 iterations, except for one stochastic orbit, which is plotted for 30,000 
iterations to improve the definition of the stochastic web. The vertical axis is the dimensionless 
angular momentum (^) and the horizontal angle is the apocenter azimuth. The dashed lines show 
the maximum dimensionless angular momentum allowed by (^). (b) The map (p)^-|T7|) with 

parameters a = 0, m = 2, e = 0.5. Resonant orbits of period 2 (“bananas”), 4 (“fish”), and 
6 (“pretzels”) are marked by circles, squares, and triangles on the map; open and filled symbols 
denote different orbits. 
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near-circular orbits are present in the integration but not the map (between the last plotted orbit 
and the dashed boundary, centered on (j) = ^tt, ^tt); (iii) in the map there is no maximum angular 
momentum, so that some orbits continue to values of |y| > 1 (these are not shown since they 
have no physical relevance). These differences are fairly minor, and mostly occur for near-circular 
orbits, for which the map was never intended to work well. Overall, the remarkable similarity of 
the two figures shows that the map correctly captures the qualitative dynamics of eccentric orbits 
in the logarithmic potential. 

The SOS contains three main types of orbit; 


• Loop orbits: these cross 0 = 0 with angular momentum \y\ > 0.5. For each orbit sgn(y) is 
fixed; that is, they circulate around the center in a fixed sense. The boundary of the loop 
orbits is defined by a separatrix that passes through the unstable fixed points at 0 = ^vr, |7r, 
y = 0. 


Box-like orbits: Orbits with zero average angular momentum, which occupy a set of resonant 
islands inside the separatrix. The orbits are box-like because the apocenter azimuth librates 
around 0 = 0 or vr and avoids 0 = |7r and |7r. These orbits were examined by Gerhard fc 


Binney (1985), Pfenniger &: de Zeeuw (1989), and christened “boxlets” by Miralda-Escude &: 


Schwarzschild (1989). Each periodic boxlet has a signature which is the sequence of the signs 


of the angular momentum at successive apocenters (-t-, — or 0). The periodic orbits at the 
centers of the major islands correspond to “banana” orbits (period 2, signature [00]), “fish” 

orbits (period 4, signature [0-1-0—]), and “pretzel” orbits (period 6, signature [0-I--I-0-]); 

the nomenclature is that of [Binney (1982) and Miralda-Escude &: Schwarzschild (1989). 


• A stochastic web, which shows up on the SOS as a connected fuzzy structure passing through 
y = 0, 0 = ^TT, Itt. Despite the complicated geometry of the web, this is a single orbit, which 
was followed for 30,000 iterations (compared to 1000 for the other orbits) to improve the 
coverage. 


The SOS also provides preliminary information on which orbits are needed to construct 
self-consistent equilibrium galaxy models. The torques described by equation (^) with m = 2, 
or the potential (24), are generated by a density distribution that is elongated along the axis 
0 = 0, vr. The SOS shows that loop orbits are most eccentric when their apocenters are at 
0 = ^TT, Itt and hence these do not support the required elongation of the density distribution. 
The box-like orbits, however, have apocenters that librate around 0 = 0 or vr and avoid ^vr and 
Ivr; therefore these are aligned with the required figure. However, detailed studies of orbits in the 
logarithmic potential are required to determine whether the density distribution is narrow enough 
to allow self-consistent equilibrium models (Richstone 1980, Pfenniger & de Zeeuw 1989, Lees & 
Schwarzschild 1992, Kuijken 1993, Schwarzschild 1993). 


Let us now consider the stochastic web in more detail. Because of the singularity in the 
logarithmic potential at r = 0, we expect that any orbit passing close to the origin will be chaotic. 
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In the context of the map, an orbit leaving the origin has zero angular momentum and azimuthal 
angle cj)-, shortly before apocenter it receives an angular momentum kick —^esinmcj) (cf. eq. 1^); 
therefore its phase-space coordinates at apocenter are 


{(p,y) = (0o,-5esinm(/>o). 


(28) 


All points in the map satisfying this constraint—and their iterates under the map—should belong 
to stochastic orbits. What is less obvious is that they all belong to the same orbit, i.e. that all 
orbits passing through the origin are part of the same stochastic web. We cannot prove this result 
but it is consistent with all of our numerical experiments with logarithmic potentials. 

The geometry of the web is particularly striking when e is small. Figure |a shows the web 
generated by 100,000 iterations of a single orbit in a map with a = 0, m = 2, e = 0.001. A very 
similar plot is generated by direct orbit integration in the analogous non-axisymmetric potential. 


The following analysis provides some insight into the behavior shown in Figure |3|a. 
V’n = 0n — ngo{a) mod(27r). Since go{a) = n for a > 0, ipn = (kn if n is even and ijjn = (t'n 
is odd. Then for m = 2 and a > 0 equations (15-17) become 


Let 

± TT if n 


y'n = yn - |esin2V’n, 

V'n+l = 1pn+g{a,yn) - 9o{a), 

yn +1 = y'n-Iesin2ipn+1- (29) 


This mapping is equivalent to motion under the Hamiltonian 

OO 

H = G{a,y) — cos 2^lJ = G{a,y) — cos 2'ijj ^ exp{27rijt). (30) 

j=-oo 

Here Si{t) = “ ^) is the periodic delta-function with unit period, t is a continuous 

variable which equals n at the iteration of the mapping, y is the momentum conjugate to the 
coordinate V', = n) = ipn, y{t) = for n < t < n -|- 1, and 

ry 

Gia,y)= [g{a,y) - go{a)]dy. (31) 

Jo 

If the motion is slow, the terms in (|30D with j ^ 0 have relatively little effect, and we can 
approximate the Hamiltonian by its averaged value 

H = G{a,y) — ^ecos2'ijj. (32) 

The motion is along level surfaces of this Hamiltonian, which are plotted in Figure ^3 for the 
same values of a and e as the map in Figure ^a. We see that the averaged Hamiltonian correctly 
captures the location of the separatrix and the division into loop orbits and box-like orbits (but 
not, of course, the detailed structure of the stochastic web). 
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- 0.002 
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Fig. 3.— (a) A single stochastic orbit in the map (p!^-|T7|) with parameters a = 0, m = 2, e = 0.001. 
The orbit has initial conditions 4> = ^vr, y = —0.0005 (cf. eq. |2^) and is iterated 100,000 times, (b) 


Level surfaces of the averaged Hamiltonian (32), for the same parameters as in (a). The variable ■0 
equals cj) or 4> ±Tr. Level surfaces corresponding to loop orbits are shown by solid lines, and those 
corresponding to box-like orbits are shown by dotted lines. 
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2.4. The Liapunov exponent 


The map can be used to determine the Liapunov exponent associated with the web for any 
value of the perturbation strength e. We simply follow the tangent map for an orbit that passes 
close to the center, i.e. one with initial coordinates given by equation (28). The tangent map is 
defined by taking differentials of equations (1^-lT): 


The Liapunov exponent is 


dy'n 

1 

e 

II 

^me cos{m(pn)d(pn, 


d(pn+l 

— d(pn “1“ 

Qy 1*^’ Vnldyni 


dyn +1 

1 

e 

II 

ime cos mcpn+id(pn+i ■ 

(33) 

A = 

: lim — In 

f d(pl + dyiy^‘^ 

(34) 


n^oo 77, 

1 d(pl + dyl j 


In practice we estimate the Liapunov exponent using large but hnite n. Since we expect that all 
initial conditions (28) he in the web (i.e. they are all the same orbit), the Liapunov exponent 


should be independent of the initial azimuth (pQ. 


Figure ^ shows the Liapunov exponent A for maps with a = 0 and m = 2. The initial azimuth 
4>q was chosen randomly in the interval [0, 27r) and each map was iterated 100,000 times. Although 
there are some outliers, in general the scatter is small, confirming that A is almost independent 
of (pQ. What is unexpected is that A depends only weakly on the perturbation strength e—from 
around 0.07 near e = 0.3 to 0.02 near e = 10““^. 


3. Concave potentials (a > 0) 

The maps with m = 2, e = 0.3 and a = 1.5,1, 0.5 are shown in Figure ^a,b,c. Like the map 
for the logarithmic potential examined already, these exhibit three types of orbits: loop orbits, 
which have a fixed sign of angular momentum; box-like orbits, which have zero average angular 
momentum and an apocenter azimuth that librates about (/> = 0 or tt; and a separatrix orbit, 
which passes through y = 0, cp = ^tt, and separates the loop orbits from the box-like orbits. 

Panel (d) shows the SOS for orbits integrated in the potential (^) with a = 1, b = 0.91, 
which corresponds to the map (b). The similarity of the map in panel (b) to the SOS in panel 
(d) conhrms that the map captures the qualitative dynamics in the non-axisymmetric potential. 
The principal difference, apart from minor changes in orbit shape, is the small resonant islands for 
near-circular orbits that appear in the integration but not the map (near the dashed boundary, 
centered on <p = ^vr, |7r). 
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e 


Fig. 4.— The Liapunov exponent A (eq. as a function of perturbation strength e, for maps 
with a = 0 and m = 2. The initial conditions were taken from ( |^ ) with cpQ chosen randomly in 
[0,27r). The perturbation strength e was also chosen at random, from a uniform distribution in 
loge between —4 and —0.5. A total of 500 {4>o,e) pairs were examined, and each was iterated 10^ 
times. 
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Fig. 5.— The map (p!^-|T7|) with parameters m = 2, e = 0.3, and (a) a = 1.5; (b) a = 1; (c) 
a = 0.5. The orbits have been iterated 1000 times, except for the orbits at the separatrices, which 
are iterated 5000 times to provide better coverage, (d) The SOS for orbits in the potential (M), 
with a = 1 and axis ratio b = 0.91 (cf. eq. 26). The dashed lines show the maximum dimensionless 
angular momentum allowed by (25). Each orbit is plotted for 1000 periods, except for the separatrix 
orbit, which is plotted for 5000 periods. 
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The locations of the separatrices are well-described by the averaged Hamiltonian (32), as 
illustrated in Figure 

There are two interesting transitions in the behaviour of the box-like orbits as a is varied: (i) 
In Figures |^,b {a > 1) most of the box-like orbits form smooth curves surrounding the points 
y = 0, (/> = 0, vr, while in Figures |^c and (0 < a < 1) the box-like orbits are almost all broken up 
into resonant islands. Miralda-Escude &: Schwarzschild (1989) call orbits of the first type “normal 
boxes” and those of the second type “boxlets”. Presumably this transition reflects the validity of 
the standard map as an approximation to the present map when a > 1 (i3). (ii) In Figure § 
the separatrix orbit is part of a stochastic web that surrounds all of the resonant islands, while in 
Figure the separatrix orbit does not penetrate between the islands. Presumably this difference 
arises because KAM surfaces divide the separatrix from the islands in the latter case. 


In §2.1. we argued that for 1 < a < 2, the mapping resembles the Taylor-Chirikov map, with 
the critical e for global stochasticity given in equation (20). Thus the motion should be nearly 
regular for small e, in agreement with the phase portraits shown in Figures |^,b and d. 


3.1. A simpler map 


The behavior of box-like orbits in concave potentials can be clarified by examining a simpler 
map. We begin with the map (^) and assume that V' is small, so that sm2^jJn — We also 
replace g{a,y) — go{a) by its asymptotic form for small y (eqs. ^ and HD- Thus we may write 

y'n = Vn- e'ljjn, 

'ijjn+i = -ipn + C sgn^y'Jly'J^, 

Vn+l = y'n-e'll^n+l- (35) 

where C > 0 is a constant and 0 < /3 < 1. We now re-scale the coordinates, choosing new variables 
/ = (eC)V(/3-i) y' and 0 = these coordinates the map takes the simpler form 


I'n = 
^n+1 = 

In+1 — 

The motion is governed by the Hamiltonian 

H = 


In Gni 

Bn + Sgn{fn)\l'nf, 

In — Gn+lj 


7 


+ 6i{t)e^ 


(36) 


(37) 


where 5i{t) is the periodic Dirac delta function with unit period, j = (5 + 1, 9n is the value of 
9{t) at t = n, and In is the value of /(t) for n < t < n -|- I. We consider first the time-averaged 
Hamiltonian 

Ho = \I\Vl + 9^ 


(38) 
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Fig. 6.— Separatrices of the averaged Hamiltonian (^), for m = 2, e = 0.3, and a = 1.5,1,0.5. 
These agree well with the locations of the separatrix orbits in the corresponding maps, panels 
(a),(b),(c) of Figure |. 
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The corresponding phase-space portrait consists of closed curves that are not differentiable on the 
axis / = 0. With each trajectory we associate the action: 


J = — / IdO 
2ttJ 


TT 


0 > 


(39) 


where — x^)^dx = 7 r^/^r(l -|- l/ 7 )/r(| -|- I/ 7 ). The energy is simply given by; 


2 ^ 


Hq{J) = cvqJ' 1 +'^, with 


1 


TT 


27 

7+2 




(40) 

77+2 \” 1 / 

Note that the frequency oj{J) = dHo/dJ diverges as J tends to zero, since 7 = /3 -|- 1 < 2. This 
feature implies that the motion is highly unstable near the origin once we add the high-frequency 
terms of the periodic delta function. 

The transformation to action-angle variables is accomplished with the generating function 


5(0,j) = 7 “ / [Ho{j)-e^]-de, 


(41) 


which naturally gives the angle conjugate to J, ip = dS/dJ. 

Now we restrict ourselves to the case /3 = 7 — 1 <C 1, which is valid for potentials that are 
close to the logarithmic potential. We have 7 ^ 1 , Ho{J) = cvq = ( 37 r)^/^ 2 “"^/^. We can 

also write 6 = {dH^/dJ)~^g{(p), where g{(l)) is the 27r-periodic sawtooth function defined by 


g{x) = <^ TT-X, 

g( 27 r + x). 


|x| < ^TT, 

^TT < X < ItT, 


(42) 


In the action-angle variables (J, (p) of the averaged Hamiltonian, the mapping Hamiltonian 
becomes 

H = HoiJ) + Mt) - 1] \\<P) = + ^Mt) - l]ff2(0)| . (43) 

We would like to analyze the principal resonances of this Hamiltonian. After Fourier 
expansion of g and 5, the Hamiltonian takes the form 


H = HoiJ) 


2 “ , , 4 ^ “ (- 1 )^ 

i+ji:cos(2™()+^ee^ 

n=l rn^0n=l 


cos 2im(p — Trnt) 


(44) 
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Choosing a resonant pair {m,n), and including all terms of the form (pm,pn), we end up with the 
resonant Hamiltonian: 


Hmn = HoiJ) 


4 “ (-1)^’ 

1 H- ^—2 2^ — Z3 .— 2p{m<p — imt) 


m^TT 


p^ 


(45) 


We change coordinates with the canonical transformation: F[J\(j)) = 2J'{m([) — nrit), giving 
J = 2m J' and cj)' = 2{m4> — Trnt), and end up with the one degree of freedom Hamiltonian 


which simplifies to 


1 + 


E 


(-i)p" 


2 2 
tZ/tt 

P=1 


■ cos{p4>') 


— 2'KnJ', 


(46) 


H'mAJ'A') = Hoi2mJ') 
= Ho{2mJ') 


1 + 


1 - 


2 <^' ^ 


3m^ m^vr m^vr^ 


— 2'KnJ', 0 <(?!)' < 27r, m even 


1 


a /2 


+ 


3m^ m^TT^ 

This Hamiltonian has equilibria at 0 and vr. For m even and 7 ~ 1 we get 

2m7a;o(l + 


— 2'KnJ', —TT < (f)' < TT, m odd. (47) 


Jo = 


Jtt - 


2 

3m^ 


2+7 

2-7 


7rn(7 + 2 ) 

2m7a;o(l — 


1 

3m^ 


_ 2+7 
\ 2-7 


7rn(7 + 2 ) 


(48) 


and for m odd, Jq and exchange values. The difference between Jq and Jjr diminishes with 
increasing m. The resonant tori in 1-9 space are recovered via equation (^). 

We observe that, for /3 -C 1, the stochastic web outlines the separatrices of the resonances 
corresponding to (m > 0, n = 1). In Figure 0, we show, in solid lines, the resonant tori 
corresponding to J,r and (m = 1... 15, n = 1), for /3 = 0 . 2 . Superposed on them is the web 
structure that obtains at this value of /3 from the map (|3^ . The tori string together the periodic 
points in chains of islands. These periodic orbits and the “boxlets” associated with them dominate 
the phase space. The separatrices, nested in tangent contact, form the network along which orbits 
diffuse in the web. 


4. Convex potentials {a < 0) 

The maps with m = 2, e = 0.3 and a = —0.25, —0.5, and —0.75 are shown in Figure ^,b,c. 
These show the usual loop and box-like orbits, as well as a chaotic zone arising from a single orbit. 
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Fig. 7.— The resonant tori that are predicted by ( |48| ) for n = l,m = 1... 15. The web 
determined by the map (^6|) is superposed on top, along with one orbit near the center of each 
of the major islands. The predicted tori pass through the periodic points in the largest chains of 
islands outlined by the web. 
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Panel (d) shows the SOS for orbits integrated in the potential (^) with a = —0.5, c = 0.05, which 
shows good qualitative agreement with the map in panel (b). Note that the stochastic web seen 
near a = 0 develops into a stochastic sea by a = —0.75. 


5. Other potentials 

5.1. m = l 

We can use the map to explore the behaviour of orbits in lopsided (m = 1) potentials. In so 
doing we restrict ourselves to the range a < 1 because dipole potentials with 1 < a < 2 cannot 
be generated by a plausible density distribution (cf. eq. |2^). The maps with m = 1, e = 0.5 and 
0 = 1, 0.5, 0, and —0.5 are shown in Figure Pa,b,c,d. These show the usual loop orbits, box-like 
orbits, and chaotic orbits. 

The maps provide preliminary estimates of which orbits are needed to construct self-consistent 
lopsided galaxies. The torques described by equation (H) are generated by a density distribution 
that is elongated along the axis 0 = 0. When a = 1 (panel [a]) the loop orbits are most elongated 
(smallest |y|) when their apocenters are near 0 = 1.8,4.6 and the box-like orbits have approximate 
m = 2 symmetry around 0 ~ 0, vr; hence neither orbit family is likely to support a lopsided density 
distribution. The situation is more promising for a < 0 (panels [c] and [d]). Here both the loop 
orbits and the chaotic sea are most elongated at 0 = 0. Thus it may be possible to construct 
self-consistent lopsided galaxies with centrally concentrated density distributions. 


5.2. Potentials that are not scale-free 

In writing the precession rate as g{a,y) (eq. |^, we have assumed that the underlying 
axisymmetric potential is scale-free. More generally the precession rate is a function of the energy 
of the orbit in the underlying axisymmetric potential, g{y,E), where E = -|- <l>(r) is assumed 

to be conserved since the non-axisymmetric component of the potential is small. Apart from the 
energy-dependent precession rate the map would remain unchanged. 


5.3. Rotating potentials 

The map can also be used to study orbits in rotating potentials. Let us suppose that the 
non-axisymmetric component of the potential is stationary in a frame rotating with pattern speed 
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Fig. 8.— The map (HMD with parameters m = 2, e = 0.3, and (a) a = —0.25; (b) a = —0.5; 
(c) a = —0.75. The orbits have been iterated 1000 times, except for the chaotic orbits, which 
are iterated 5000 times to provide better coverage, (d) SOS for orbits in the potential (^), with 
a = —0.5 and c = 0.05. Each orbit is plotted for 1000 periods, except for the chaotic orbit, which 
is plotted for 10000 periods. 
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Fig. 9.— The map (MD with parameters m = 1, e = 0.5, and (a) a = 1; (b) a = 0.5; (c) 
a = 0, (d) a = —0.5. The orbits have been iterated 1000 times, except for the chaotic orbits, which 
are iterated up to 20000 times to provide better coverage. The range 1 < a < 2 is not examined 
because dipole potentials with this radial exponent cannot be generated by a scale-free density 
distribution (cf. eq. |2^). 
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Qp, that the potential has m-fold symmetry and that its azimuthal minimum lies along the ray 
(p = Opt. The azimuthal angle in the rotating frame may be written ip = p — Qpt. Then if 
Tj.{y,E) is the radial period of the particle orbit in the underlying axisymmetric potential, the 
map (H) -(p^ is easily generalized to 


y'n = Vn- sin mijjn, 

V^n+l = 'Pn ~i~ siVn, E) — ^pTr{yji, E), 

Un+I = y'n - sin mipn+i- (49) 


Note that g is an odd function of y while is even. 


5.4. Three dimensions 


We derive a mapping for motion in the three-dimensional generalization of the potential 


sgn(a) , a 7 ^ 0, 

\ log + -^ + ) 0 = 0. 


(50) 


Here an orbit is defined by its vector angular momentum L and the unit direction vector 
of its apocenter, n, where L • n = 0. We again prefer to deal with the dimensionless angular 
momentum, so we define Y = L/Lc(E). The torque gives a kick to the angular momentum vector 
which is concentrated at apocenter and equal to —Qaiin x t„ where t = {nx-,ny/}p‘^Uz/c^)- For 
nearly spherical potentials, Qa is specified in (^). This kick is followed by a rotation of n about 
the updated Y. Putting it all together. 


Y' = 
rin+l — 
Y„+i = 


Yji 2^0^*^ ^ 

1) £ 

Yjj -^Qa^n+l X tn+li 


a(a lY' n 

n{Y'^,g{a, |Y;|))n„ = x Yi 


lY' I 

I n I 


(51) 


where TZ{v, 9) is a rotation about axis v by an angle 9. 

The mapping is apparently six-dimensional. However, it preserves the unit magnitude of n 
and the orthogonality of Y and n. These two constraints on the motion reduce the dimension to 
four. Also, the mapping is symplectic. The proof is deferred to Appendix B. 


A useful approximate description of the motion, valid when the potential is not far from 
spherical, is obtained by replacing the mapping equations (51) by the analogous differential 
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equations 


dY 

dn 

dn 

dn 


-Qo-n X t, 


g{a, |Y|) 

|Y| 


n X Y, 


where n is the orbit number. These have the integral of motion 


(cf. eq. [BTD 


F(Y,n) 


ry 

/ 9{a,y)dy + 2Qan-t. 


(52) 


(53) 


where y = \Y\ and H is an averaged Hamiltonian analogons to eqnation ([^). This integral can 
be used to solve for y in terms of n, so the equations of motion ( |5^ can be replaced by equations 
for the two perpendicular unit vectors n and y = Y/y: 


Wet r. ^ ^ /J. M 

= - t X n - yy (t X n)] , 

dn y 

^ = -5'(a,y)nxy. (54) 

The solutions of these equations of motion can be investigated using a snrface of section (for 
example, defined by n • 62 = 0 ), but we shall not do so in this paper. 


6. Discussion 


Onr map provides a heuristic tool for studying the behaviour of eccentric orbits in 
non-spherical potentials. The map reproduces most of the qualitative features of orbits in 
non-axisymmetric potentials similar to those found in realistic galaxy models (e.g. Fig. ^). The 
map can be applied to a variety of potentials, as outlined in §|^ In a limited sense, the map is 
more general than orbit integrations, since it does not depend on the specific radial form of the 
non-axisymmetric potential so long as the torque is concentrated near apocenter. The map is 
faster than direct numerical integration by 2-3 orders of magnitude—even more for near-radial 
orbits in centrally concentrated density distributions, which are difficult to integrate acenrately. 
The speedup offered by the map is particularly important when exploring the long-term evolntion 
of orbits in the central regions of galaxies: at 10 pc from the center of a typical giant galaxy the 
crossing time is only 10“® of the Hubble time. 

The map offers a powerful approach to studying many aspects of the behaviour of orbits 
in triaxial potentials, including the distribution of Liapunov times in chaotic orbits, the rate 
of mixing of a non-random distribution of stars, trapping of chaotic orbits, and the long-term 
inflnence of a central black hole on centrophilic orbits (Merritt 1996, Merritt & Fridman 1996, 
Merritt & Vallnri 1996). 
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It may also be instructive to use the map to study the equilibrium and stability of approximate 
stellar systems, in which the torque between two elongated orbits is determined by the relative 
orientation of their apocenters. A very similar approximation has been used as the basis of 
heuristic explanations of bar formation in disk galaxies (Lynden-Bell 1979) and the radial-orbit 
instability in spherical galaxies (Palmer &: Papaloizou 1987). 


We thank Jeremy Goodman for thoughtful advice. This research was supported in part by 
NSERC. JT acknowledges the support of the Harlan Smith Fellowship under NASA grant NAGW 
1477, and ST acknowledges the support of an Imasco Fellowship. 


A. Power series for the precession rate in power-law potentials 


The goal of this Appendix is to provide a power-series expansion for the precession rate g{a, y) 
(eq. 0 )- For the sake of brevity we focus first on potentials with a > 0 (cf. eq. ||), extending the 
results to a < 0 at the end. 


We first extend the domain of g by defining g{a,y) = 0 for y > 1. We then take the Mellin 
transform 

POO 

G{a,s)= g{a,y)y''~^dy] (Al) 

Jo 

upon evaluating the integrals we have 


^i/22./2r(i + i,)r(,/a) 

G[a, s) = N I 1 —^—TT’ " > 

which holds in the right half-plane Re(5) > 0. 

The inverse Mellin transform is 

P(7-\-ioO 


(A2) 


P(7-\-lOC 

9 {a,y) = — G{a,s)y~ 

2,7ri Ja—ioo 


Vs, 


(A3) 


where a is real and positive. For 0 < y < 1, the integration contour can be closed by a semi-circle 
in the negative half-plane from a + ioo to u — ioo. Since \/T{z) is entire, the only contributions 
to the contour integral (A3) arise from poles in the functions r( 2 ;) in the numerators of equations 
(A2), which occur at 2 ; = —n, n = 0,1, 2,... and have residues (— l)"'/r(n -|- 1). Thus 


9{a-,y) 


^ 1/2 ^ (-l)’"h^^+^(a)r[-(2n-H)/«] on+i 

2’"-V2ar[i - re - (2n + 1)/a]r(re + 1) ^ 

1/2 (-l)-h-(«)r[(l-«n)/2] 

2 “"'/ 2 r(i — n — an/2)T{n -|- 1) 


(A4) 
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Using the identity r( 2 ;)r(l — z) = Tr/sinvrz this expression can be rewritten as 


^^"■^na)cot[7r(2n + l)/a]r[|+n + (2n + l)/a] 2n+i 

— TT TT 2_^ nn — ^/0 -nr -1 , /r. , i \ / 1 -n/ , i\ V 


n=0 

^ h°‘'^{a) tan(i7ran)r[n(l + a/2)] 


V2ar[l + (2n + l)/a]r(n + 1) 

a > 0. 


n=l 


2“”/2r[(l + a!n)/2]r(n + 1) 


(A5) 


In the limit 7/ —> 0 we recover equations and (0)- 

For a < 0 we can apply the relation ( p^) to obtain 

. X ^ 27r 1/2 /t2n+i[2/?/(2 - (5)] tan[7r(2n + l)//?]r[(2n + 1)//?] 2 n+i 

^ 2-/1^ 2’^-V2/3r[i - n + (2n + l)//3]r(n + 1) ^ 

1/2 /t2/3n/(2-/3)p/j/(2 _ /?)] tan[7r/3n/(2 - /?)]r[2n/(2 - (3)] 20 n/i 2 -g) / aci 

2/3«/(2-/3)-i(2-/3)r[i + /3n/(2-/3)]r(n + l) ^ ^ ^ 


where /3 = —a, a < 0. The series ( |A5D and ( [A6P are asymptotic: if truncated after a fixed number 
of terms the series become arbitrarily accurate as y —!■ 0, but for fixed y the series do not generally 
converge as more terms are added. 


B. The three-dimensional mapping is symplectic 


We claim that the generalization of the two-dimensional mapping (eqs. |T5-17) to triaxial 
potentials (eqs. is symplectic. There are two routes to the proof: one familiar but messy, 
the other less familiar but more elegant and illuminating. We outline the first and describe the 
second in detail. On the messy road, we define two pairs of canonically conjugate variables. One 
pair joins the magnitude of the angular momentum and the argument of apocentre and the other 
the X 3 -component of the angular momentum vector and the longitude of the ascending node on 
the X 1 -X 2 plane. That these coordinates are canonically conjugate is a basic result of classical 
celestial mechanics (e.g. Plummer 1960). Then we express the map in terms of these coordinates, 
derive the Jacobian M of the mapping transformation, and show after some algebraic labor, 
that this Jacobian fulfills the requirements that a canonical transformation must satisfy, namely: 
M^JM = J, where J is the 4x4 canonical symplectic matrix (see Goldstein 1980 for definitions 
and details). 


We prefer to observe that the mapping derives from the Hamiltonian: 

H = Hy{Y)+di{t)Hnin), 

where Hy = p g{a,y)dy, y = |Y|, and Hn = • tj with the help of the bracket: 


(B7) 


{Fi(Y,n),F2(Y,n)} = Y • (Vy^i x VyT 2 ) + n • (V^Fi x VyT 2 + Vy^i x (B8) 
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The bracket 0 gives the time evolution of any real valued function Fi of these two vectors, along 
the vector field generated by another such function F 2 , through the relation 


|Ti(Y,n) = {Fi,F2}. 


(B9) 


The equations for the time evolution of Y and n along the vector field generated by Hy are: 

d. 


dt 

!L 

dt 


Y = 0, 


-n = —- 


5(a, |Y|) 


n X Y. 


(BIO) 


The first of these equations states that Y is conserved by Fly. The second is telling us Hy causes a 
rotation of n by angle g about vector Y. This is one of the mapping steps. Hn on the other hand 
generates motion: 

-Qall X t, 

0. (Bll) 

When we compose the action of Hy and H^, in the manner dictated by the periodic delta 
function, we recover the mapping (^Tj). Now we still have to show that the dynamics generated 
by this bracket is canonical. The bracket is known as a Lie-Poisson bracket and is closely related 
to the bracket used in the study of conservative rigid body dynamics (see Marsden 1992 for a 
general discussion, and Touma and Wisdom 1994 for a discussion of Lie-Poisson brackets in an 
astronomical context). On a suitably defined pair of canonical variables, the bracket reduces to 
the canonical Poisson bracket and the dynamics generated by the bracket consists of canonical 
transformations. Thus, Hy and H^ separately generate symplectic mappings and their composition 
via the mapping is therefore symplectic, completing our proof. 



d 
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